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ABSTRACT 


A numerical procedure is derived for the solution of the Vlasov-Poisson system 
of equations in two phase-space variables. Derivatives with respect to the phase- 
space variables are approximated by a weighted sum of the values of the distribution 
function at properly chosen neighboring points. The resulting set of ordinary differ- 
ential quatlons is then solved by using an appropriate time-integration scheme. The 
accuracy of the proposed method is tested with some simple model problems. The 
results for the free-streaming case, linear Landau damping, and nonlinear Landau 
damping are investigated and compared with those of the splitting scheme. The 
proposed method is found to be very accurate and efficient. 


1 . INTRODUCTION 


It is well known that the Vla.sov equation adequately describes the nonlinear 
evolution of collisionless plasmas. Since knowledge of the nonlinear behavior of 
plasmas in two and three dimensions is indispensable in understanding the plasma 
physics of controlled thermonuclear fusion, the numerical integration of the Vlasov 
equation has been studied Intensely during recent years (refs. 1-13). However, little 
proj'.ress has been made on the development of fast and accurate integration schemes 
for the Vlasov equation in two and three dimensions for a magnetized plasma. A spllt- 
tlnu scheme by Cheng and Knorr (ref. 12), in which the Vlasov equation Is integrated 
In the original phase space by splitting the convective and acceleration terms in such 
a way that the overall scheme is second-order accurate in At, has been successfully 
applied and is one of the most promising schemes today. In two and three dimensions, 
however, the interpolation methods used in the scheme become more and more complicated 
and time-consuming, especially for magnetized plasmas (refs. 13 and 14). Development 
of new, accurate and efficient methods is still needed before the Vlasov equation can 
be solved for three-dimensional magnetized plasmas in reasonable computation time. 

In this ,'>aper a new numerical method called the Modified Differential Quadrature 
(M.D.Q.) method, is proposed to integrate the Vlasov equation. The new method is an 
extension of the Differential Quadrature (D.Q.) method proposed by Bellman et al. 

(ref. 15). In the present method, derivatives of the distribution function with 
respect to the phase-space variables are approximated by a weighted sum of the values 
of the distribution function at properly chosen neighboring points to generate a set 
of ordinary differential equations in time, whereas in the original D.Q. method, 
these are approximated by using values at all mesh points in the computat ional domain. 
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As a rssulti coaputstionsl sfflclsncy Is slgalficsntly Isprovsd with ths N.D.Q. 
method. By changing the weighting coefficients, the spatial derivatives can easily 
be approximated with as much accuracy as desired. The resulting set of ordinary dif- 
ferential equations is then Integrated by using an appropriate time- Integration 
scheme. This solution process gives a very accurate and flexible method which Is 
simple and straightforward to program. The present method has the feature that It Is 
accurate to arbitrary order In space by changing the freighting coefficients and " 

In time by choosing a suitable time- Integration scheme. Although In this paper, che 
method Is presented only In one dimension In order to Illustrate Its basic elements 
clearly. Its extension to two and three dimensions Is straightforward. 

Section 2 describes the M.D.Q. method for the Vlasov equation and Section 3 
describes the procedure for determining the weighting coefficients. In Section 4 we 
demonstrate the accuracy and efficiency of the present method through numerical 
experiments on simple model problons. In Section S we present the numerical results 
obtained for the free-streaming case, linear Landau damping, and nonlinear Landau 
damping. These results are then compared with those obtained by the splitting scheme. 
Section 6 presents the conclusions. 


2. MODIFIED DIFFERENTIAL QUADRATURE METHOD 


The system of equations under consideration consists of the one-dlmenslonal 
Vlasov equation for the electron distribution function f(x,v,t) 
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3t 




( 1 ) 


and Poisson's equation for the electric field E 


3E 

3x 



( 2 ) 


These equations are written in dimensionless form. The basic units of time t 
and velocity v are the reciprocals of the so-called electron plasma frequency u 
and the electron mean thermal velocity v^, respectively. The dis*:ance x is meas- 
ured in units of the Debye length Xy. A periodic boundary condition in x is 
assumed. A rectangular mesh will be used to represent the x-v phase space with the 

i V_.j.), giving the set of points 


(x,v) |0 i x < L, 
is the spatial periodic length and V„ax the cutoff 


computational domain R ' 

(lAx, jAv) 6 R where L 
velocity, while 1 ■ N and J - 2M designate the number of mesh points used along the 
coordinates x and v, respectively. 

It the distribution function f satisfying Eq. (1) is sufficiently smooth, we 
can write the approximate relations 
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where Index 1 stands for and J for Vj. In this paper we have modified the 

approximate relations, Eqs. (3a) and (3b), to use the values of f at the nearest 
Np and Mp mesh points centered around m and vj, respectively. Instead of using 
those at ail mesh points In the computational domain, as Is the case In the original 
D.Q. method (ref. 15). By doing this, the number of arithmetic operations to be per- 
formed for every point Is significantly reduced and, moreover. In the case of a uni- 
form mesh, the weighting coefficients aj^)^ and ajk become Independent of the Indices 
1 and j. Therefore, the approximate relations, Eqs. (3a) and (3b), can be written as 
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where a,^ - a - (Np + l)/2 In Eq. (4a) and a^ - ag^. 6 - (Mp + l)/2 In 

Eq. (4b), respectively. 


Substitution of Eqs. (4a) and (4b) into Eq. (1) yields the set of Nx2M ordinary 
differential equations in time 




dt 


^ ■ f/xi,Vj , ^ + k - a,j), ^ aj^f(i,j + 

\ KM K«l 


k - P), E(f), t 
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(5) 


The nunorlcnl solution of such a system, Eq. (5), Is a simple task using a standard 
scheme for ordinary differential equations. In the present paper, we tested three 
schemes, namely, leapfrog, corrected leapfrog, and Stetter's method (ref. 16). All of 
these schemes are second-order accurate in time. As applied to Eq. (5), tne schemes 
take the following forms: 


1) leapfrog 


,n+l 


f""‘ + 261?” 
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2) corrected leapfrog 


^n+l _ ^n-1 ^ 
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£( 0 ) , ^ 



- f" + At(F^°^ + f")/ 2 
£<2) . £« ^ At(F^^^ + f”)/2 
£"■*■* m + (i-e)f^^^ , 0 < e < 1 


where the superscripts h and n+1 refer to tine t and t+At, respectively, and the 
subscripts 1 and j are omitted. The electric field E Is computed by a standard 
direct Poisson solver technique from f by Eq. (2). 


3. DETERMINATION OF WEIGHTING COEFFICIENTS 


In order to appropriately determine the weighting coefficients in the approximate 
relations, Eqs. (3a) and (3b), the test function is taken to be the following form, by 
analogy with Lagrange's interpolation formula, 


Pk(x) * P(x)/t(x - Xj^)P'(Xj^)] 


(9) 


where P(x) is a polynomial of degree N 


P(x) = (X - Xj)(x - X 2 ) . . . (X - Xjj) (10) 

It follows that P|^(x) is a polynomial of degree N-1 such that Pjq(x£) * and 

P(xj) = 0, If the values of a function g(x) are known at N mesh points 
X =^xi, Xj, . • • XN, the polynomial of degree N-1, g(x), which coincides with g(x) 
at these collocation points can be written as 


g(x) 


g P^(x)g(x^) 


By differentiating Eq. (11) with respect to x, we have the relation 


( 11 ) 


g’ (x) - p^(x)g(xjj) 


( 12 ) 


Using the fact that such a relation as Eq. (3a) is to be exact for f(x) « Pj^(x), we 
see that 

*ik * P'(x£)/[(Xj^ - Xj^)P'(x^^)l , i k (13) 
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For the case 1 ■ k, use of 1* Hospital's rule gives 


*kk" 


In their D.Q. method, Bellman at al. have chosen as P(x) the shifted Legendre 
polynomial of degree N, pJCx), selecting x^ to be the roots of P^Cx) and 
determining a^^j^ explicitly (ref. IS). In this paper, as described In the previous 
section, we have modified the approximate relations to make use of appropriately cho- 
sen neighboring points instead of using those at all mesh points, and have adopted a 
uniform mesh spacing to make the weighting coefficients Independent jf Index 1. 


The weighting coefficients 
nlng of the calculation. 


are computed and stored once and for all at the begln- 


4. NUMERICAL RESULTS FOR MODEL PROBLEMS 


Before application to the Vlasov equation, the present method was applied to some 
simple model problems to demonstrate Its accuracy and efficiency. First we checked 
the effects of the order of approximation to spatial derivatives on the accuracy of 
computed results. Using the test function of the form 


1 + A cos mx 


In which A ■ m ■ 0.5, we calculated 3f/.'x by the approximate relation, Eq. (4a), 
for several values of Np. The spatial computational domain was of periodic length 
I, - 4 ji, and the total number of mesh points N was 16. The results are displayed in 
Table 1 in which the error is defined to be 


ERROR 


it -‘hi 


where Is the difference between the analytic solution and the numerical solution 

at the mesh point Xj^. In Table 1 we can see that as Np increases, the error 
decreases significantly, as expected. In particular, the error for Np • 7 Is at 
least 3 orders of magnitude smaller than for Np ■ 3, while the difference between 
Np ■ 11 and Np • 15 Is negligible due to the accumulation of roundoff errors in 
32-bit, single-precision arithmetic. 

The present method was then applied to the following simple linear model equation 

+ c — • 0 (17) 

3t 3x 

with c • n/4. The initial condition was 


In Table 1 we can see that as 


increases, the error 
rror for N_ • 7 is at 


f(x,0) 


1 -f A cos mx 


5 


ORttslUAL PAGE IS 
OF POOR QUALITY 

with A ■ m - 0.5. The enelytlc solution of Eq. (17) Is known to be 


£(x,t) ■ 1 + A cos n(x - ct) 


( 19 ) 


irtilch gives e meens'of checking the sccurecy of the nuaerlcsl results. The compute- 

tionsl conditions are the same as for Eq. (15). 

% 

Results are presented^fur^two aeasures of the accuracy. One la F(0.5), the 
fundauental cosine node of f Ilk tine t, defined as 
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F(0.5) ■ J ^ cos mx dx 

'o 


and the other Is the L, error of f defined by Eq. (16). Both are shown In 
Figs. 1-3 for the period from t • 32 to t ■ 48. The analytical solution for F(0.5) 
is shown as a solid line, and the numerical solutions by symbols. The L 2 errors 
of f from the numerical solutions are also shown by symbols with the scale along 
the right ordinate. 

Figure 1 compares the results obtained with Np ■ 15 for three different time- 
integration schemes described in Section 2. Although the three schemes produce vir- 
tually the same results for F(0.5), Stetter's scheme gives the most accurate result 
for L 2 error. Therefore, we decided to adopt Stetter's scheme as the standard 
time- integration scheme hereafter in this paper. The results with Np ■ 3, 7, 11, 
and 15 are compared In Fig. 2, which shows the effects of order of approximation to 

the spatial derivatives. Figure 2(a) shows the results for Np ■ 7, 11, and 15; no 

significant difference can be seen for F(0.5), while L 2 error tends to decrease as 
Np increases. The difference between Np ■ 11 and Np ■ 15 is again negligible due 

to the roundoff errors. Figure 2(b) shows results for the case of Np ■ 3, which 

corresponds to the usual second-order, centered, finite-difference approximation to 
3f/3x. In this case the result for F(0.5) shows a large phase error from the ana- 
lytic solution, which shows the need for using a higher-order approximation for spa- 
tiai derivatives. The results of the present method with Np ■ 15, Stetter's time- 
integration scheme, and those of the splitting scheme of Cheng and Knorr (ref. 12) 
are shovm in Figs. 3(a) and 3(b). Although the curves for F(0.5) obtained by the 
two methods agree very well with the analytic solution, the L 2 error of the present 
method is much smaller than that of the splitting scheme. 


5. RESULTS FOR THE VLASOV EQUATION 


In this section we will present the results of integrating the one-dimensional 
Vlasov equation and demonstrate the accurccy and efficiency of the present method. 

As described in Section 2, we will use a rectangular mesh to represent the x-v phase 
space with the computational domain R ■ {(x,v)|0 < x < L, |v| s V^^x^* Throughout 
the following examplea, the cutoff velocity, V^^* taken to be 5.0. The symbols 
N and 2M designate the number of mesh points used along the x and v coordinates, 
respectively, while Np and Mp denote the degree of M.D.Q. in each direction, cor- 
responding to Np - 1 and Mp - 1 order of accuracy, respectively. 


6 


In the flrec example %rc show results for the free-streaming case with E(x.t) ■ 0 
In Eq. (1). The Initial condition Is 


f(x,v,0) ■ fo(v)(l + A cos mx) 


( 20 ) 


with A ■ m ■ 0.5 and fo(v) ■ (2it) ^exp(-v^/2). The analytic solution corresponding 
to this Initial condition Is given by 


f(x,v,t) ■ fo(v) ^1 + A cos m(x - vt)J (21) 

Here we define the density as 

p(x,t) ■ f f(x,v,t)dv (22) 

Equation (21) combined with Eq. (22) leads to the following expression for the 
density 


o(x,t) ■ 1 + exp(-m^t^/2)A cos mx 


(23) 


Numerical calculation of Eq. (22), on the other hand, gives the following expression 
for the density 


p(x,t) • tv fo(Vj) 1^1 + A cos m(x-Vjt)J (24) 

The numerical results are shown in Fig. 4 where the computed densities at x ■ 0 
(curve A) and x ■ L/8 (curve B) are plotted against time. As can be predicted from 
Eq. (23), p 1 as t and the curves in Fig. 4 tend asymptotically to 1. The 

right-hand side of Eq. (24) is the sum of periodic functions of time, which results 
in a quasi-perlodic behavior for p(x,t) called recurrence. In the present calcula- 
tion with A and M ■ 20, the predicted recurrence time 

is Tg * 2n/(mAv) > 49.00. The numerical results presented in Fig. 4 agree very well 
with this value. 

The second example shown in Fig. 5 tests linear Landau damping for the same 
initial condition as Eq. (20), but in this case, the electric field E is retained 
in Eq. (1). We used A* 0.01, m** 0.5, N ■ 8, M* 20, and tt « 1/8. The abscissa 
is the time nondimensionallzed by a)p~^ and the ordinate shows the first Fourier 
mode E(0.5) of the electric field E. The solid curve in Fig. 5 has been obtained 
by the present method with Np ■ 7, Mp * 7, and Stetter's scheme for time integration, 
while the dashed curve has been obtained by the splikting scheme. In this case the 
electric field decays exponentially and agreement of the numerljal results with 
Landau's theory is quite good up to t • 40 • except for t, slight deviation in 

Che dashed curve shortly after t ■ 32 Up~*> The recurrence effect occurs at 
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t • 46.63 In the present method and at t > 46.50 Up~^ In the splitting 

schesM, and the times are comparable to the theoretical value • 49.00 obtained 
for the free-streaming case. 

The third example tests nonlinear Landau damping. The effect of a strong non- 
linear perturbation, A • 0.5, is shown in Fig. 6 for m - 0.5, in which the evolution 
of the first three Fc<urler modes of the electric field E is shown. This problem 
has been solved by many authors (refs. 6, 10, and 12) and, therefore, is appropriate 
for evaluating th«% present method. We used N 16, M ■ 64, and At ■ 1/8. The solid 
curve and the dashed curve in the figure show the results of the present method and 
the splitting method, respectively. In this case we used the present method with 
Np “ 15 and Mp ~ 7. Initially, the first mode damps much more than predicted by the 
linear theory while the second and third modes damp much less than Landau's theory. 
After t ■ 15 <ill modes grow exponentially until t ■ 40 where satura- 

tion occurs. The results of the first and second Fourier modes obtained by the pres- 
ent method are both qualitatively and quantitatively equal to the results obtained by 
the splitting scheme, except for a slight difference after t - 50 For the 

third mode, the results show considerable differences between t ■ 10 and 

t ■ 30 «p" . We can see, however, that the overall behavior of the t‘u'*‘ee Fourier 
modes by the present method is still the same as that of the splitting arheme. 

With the mesh points used in this case, total time of execution (CPU time) was 
55 sec using a Fujitsu FACOM M-200 computer. In order to demonstrate the accuracy 
and efficiency of the present method, we compare the results obtained by the presen** 
method with those of the splitting scheme with twice as many mesh points in each 
direction. In the former method we used N ■ 16 and M ■ 64, while in the latter, 

N ■ 32 and M ■ 128. The initial conditions as well as the other parameters are the 
same as in Fig. 6. The time evolution of the first Fourier mode of the electric 
field, E(0.5), is plotted in Fig. 7, in which the solid curve shows the result of the 
present method and the dashed curve shows that of the splitting scheme. Although the 
difference between the two curves at times later than t ■ 50 Up~' is slightly 
larger than that noted in Fig. 6, the agreement can be Judged as excellent. Actually 
the difference in Fig. 7 is much less than that existing between the splitting method 
itself with N - 16, M - 64, and N ■ 32, M ■ 128. We can see from Fig. 7 that with 
the present method, results with nearly the same accuracy can be obtained using half 
the number of mesh points along each direction. Although the present method requires 
slightly more computation time per mesh point per time step, we can conclude that the 
present method is more efficient than the splitting scheme even for the one- 
dimensional case. For multi-dimensional cases, the present method requires less 
computation time than the splitting scheme with the same number of mesh points. 


6. CONCLUSIONS 


In the present paper we have developed a numerical procedure, based on the 
modified differential quadrature method, for the solution of the one-dimensional 
Vlasov-Poisson system of equations. The derivatives with respect to the phase-space 
variables in the Vlasov equation are approximated as a weighted sum of the distribu- 
tion function at properly chosen neighboring points. The weighting coefficients are 
determined similarly to Lagrangian interpolation. As a result, the approximation to 
phase-space derivatives can be of arbitrary order of accuracy by changing the number 
of mesh points used in the approximating relations. The time integration of the dis- 
tribution function is performed as in the case of ordinary differential equations. 
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We have ehown eome exaaplea chat demonttraCt the accuracy of the preaant 
numerical method. The reaultt for atrong nonlinear Landau damping are In excellent 
agreement with Choie ualng the splitting scheme. Using half the number of mesh points 
In each direction, we have found that cha present method is as accurate as the split* 
ting method. The Important feature Is that tha computation time is not only quite 
luw, but also that accurate results can be obtained with a fewer number of mesh points 
due to the high accuracy of approximation to the spatial derivatives. 

Although we have shorn the accuracy and efficiency of Che present method only In 
one-dlmenslonal problems, the extension to the two* and three-dimensional Vlasov equa- 
tion is simple and straightforward. Compared with the splitting scheme, the present 
method Is more efficient In multl-dlmenslonal problems. In fact, for the two* 
dimensional Vlasov equation, we have obtained accurate solutions In less computation 
time than the splitting scheme with the saiee number of mesh points. The results will 
be presented in a forthcoming paper. 
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TABLE 1 L 2 errors OF THE M.D.Q. METHOD WITH Np • 3, 7, 11, and 15 for 

FUNCTION EQ. (IS) 


9/ ax OF THE 


Np 

3 

7 

11 

15 

Error 

0.45087E-2 

0.46223E-S 

0.55815E-6 

0.52457E-6 
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Fig. 1 Comparison of accuracy between three time 
integration schemes for a linear model equation. 
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Fig. 2 Comparison of accuracy between the M.D.Q. 
methods with different orders of spatial approx- 
imation. (a) Np ■ 7, 11, and 15; (b) Np • 3. 
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Fig. 4 Time evolution of the density at positions 
X ■ 0 (curve A) and x - L/8 (curve B) , for the free- 
streaming case, and for the initial condition in 
Eq. (20). A - m - 0.5, N - 8, M - 20, Np - 7, and 
- 1 / 8 . 
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Fig. 5 Electric field vs. time for linear Landau damping 
with A - 0.01, m - 0.5, N ■ 8, M • 20, Np - 7, Mp • 7, and 
At - 1/8. 
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